//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

#define UNCORRELATED

#include "common.cpp"

static void usage (const char *argv0)
{
  fprintf (stderr, "%s: %s [-e exponent] [-n] [-r seed] [-R rnd.bin] [-o freq.csv] [-s samplecount] [-i start end] [-g generator] [-b] [-l growthpow startlayer endlayer]\n", argv0, argv0);
  fprintf (stderr, "    -e exponent - set the correlation exponent\n");
  fprintf (stderr, "    -n - non-progressive calculation\n");
  fprintf (stderr, "    -r seed - set the random seed\n");
  fprintf (stderr, "    -R rnd.bin - random binary bitstream data\n");
  fprintf (stderr, "    -o freq - output data to file freq.N.csv\n");
  fprintf (stderr, "    -s samplecount - how many samples to do at first level\n");
  fprintf (stderr, "    -i start end - calculate with A from start to end\n");
  fprintf (stderr, "    -g generator - use the given random generator, put 'list' to see available\n");
  fprintf (stderr, "    -b - add baseline\n");
  fprintf (stderr, "    -l growthpow startlayer endlayer - specify the layer growth base and layer range\n");
  exit(-1);
}

int main (int argc, char *argv[])
{
  big_t SAMPLES = 16384;
  bool progressive = true;
  const char *outfn = "freq";
  int start = 32;
  int end = 2048;
  int rnd_seed = 0;
  EXPONENT = 0.75;
  bool add_baseline = false;
  const char *rndpath = 0;

  const gsl_rng_type * rng_type = gsl_rng_default;

  for ( int i = 1; i < argc; i++ )
    {
      if ( argv[i][0] != '-' || argv[i][1] == 0 || argv[i][2] != 0 )
          usage (argv[0]);

      switch (argv[i][1])
        {
          case 'i': {
            if ( i + 2 >= argc )
                usage(argv[0]);
            if ( sscanf (argv[++i], "%d", &start) != 1 ||
                sscanf (argv[++i], "%d", &end) != 1 ||
                start < 1 || end < start )
                usage (argv[0]);
            break;
          }
          case 'n': {
            progressive = false;
            break;
          }
          case 'b': {
            add_baseline = true;
            break;
          }
          case 'o': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            outfn = argv[i];
            break;
          }
          case 'r': {
            i++;
            if ( i == argc || sscanf (argv[i], "%d", &rnd_seed) != 1 )
                usage (argv[0]);
            break;
          }
          case 'R': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rndpath = argv[i];
            break;
          }
          case 'e': {
            i++;
            if ( i == argc || sscanf (argv[i], "%lf", &EXPONENT) != 1 )
                usage (argv[0]);
            break;
          }
          case 's': {
            i++;
            if ( i == argc || sscanf (argv[i], "%llu", &SAMPLES) != 1 )
                usage (argv[0]);
            break;
          }
          case 'g': {
            i++;
            if ( i == argc )
                usage (argv[0]);
            rng_type = 0;
            for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
              {
                if ( strcmp ((*t)->name, argv[i]) == 0 )
                  {
                    rng_type = *t;
                    break;
                  }
              }
            if ( !rng_type )
              {
                printf ("Random generators available:\n");
                for ( const gsl_rng_type **t = gsl_rng_types_setup(); *t != 0; t++ )
                    printf ("  %s\n", (*t)->name);
                exit(0);
              }
            break;
          }
          case 'l': {
            if ( i++ == argc || sscanf (argv[i], "%lf", &BASEPOW) != 1 ||
                i++ == argc || sscanf (argv[i], "%d", &LAYERSTART) != 1 ||
                i++ == argc || sscanf (argv[i], "%d", &LAYEREND) != 1 )
                usage (argv[0]);
            break;
          }
          default: {
            usage (argv[0]);
            break;
          }
        }
    }

  if ( 2 * end > FIELD_SIZE )
    {
      fprintf (stderr, "Error: interval too large for the compile time field size!\n");
      abort();
    }

  rndf.seed (rng_type, rnd_seed, rndpath);

#ifdef UNCORRELATED
  printf ("Uncorrelated, ignoring the correlation exponent\n");
#else
  printf ("Correlation exponent is %f\n", EXPONENT);
  printf ("Cumulation %d\n", CUMULATION);
  printf ("Layer size growth power %g\n", (double)BASEPOW);
  printf ("Adding baseline: %s\n", add_baseline ? "YES" : "NO");
#endif
  printf ("Field size %d\n", FIELD_SIZE);

  Field *f = new Field(add_baseline);

  printf ("Running, %s from %llu samples\n", progressive ? "progressive" : "non-progressive", SAMPLES);

  for ( int A = start; A <= end; A *= 2, SAMPLES /= (progressive ? 2 : 1) )
    {
      printf ("Calculating A=%d, %llu samples..\n", A, SAMPLES);
      freq fa;

      size_t samples = 0;
      double cum_len = 0;
      while ( SAMPLES > samples )
        {
          f->reset();
          int x = FIELD_SIZE / 2;
          int y = FIELD_SIZE / 2;
#ifndef UNCORRELATED
          f->set_sealevel (f->value (x, y));
#else
          if ( !f->valueb (x, y) )
              continue;
#endif

          coord_t start;
          start.pos.x = x;
          start.pos.y = y;
          start.low = rndf.rnd() >= 0;
          Walker walker(f, start);
          if ( walker.walk (A) )
            {
              cum_len += walker.get_len();
              samples++;
              fa.add (walker.get_hull());
            }
        }
      assert ( samples == SAMPLES);

      fa.print (outfn, A);
    }
  delete f;

  return 0;
}

